R Code for GCB-19-0810

rm(list = ls())  # Remove all objects from memory

1. Influence of most intense ENSO events on the coral cover in the ETP

library(ggplot2)
library(gridExtra)

ENSO_colours<-c("blue", "red") # Type of anomaly
my_colours<-c("#feb24c93", "#fb6a4a93", "#9e9ac893", "#a6d96a93") # Thermal regime

theme_set (theme_classic() + theme(panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(), 
                              axis.line = element_line(colour = "black"),
                              legend.position="none",
                              axis.text.x = element_text(angle = 90, vjust = 0.5),
                              plot.title = element_text(size=12, face="bold"),
                              #panel.border = element_rect(colour = "black", fill=NA, size=1)
                              panel.border = element_blank()
                              ))

ENSO intensity in region 3

Figure 2a

# Plot ENOS 1971 - 2014
  ENOS<-read.table(file="http://www.cpc.ncep.noaa.gov/data/indices/ersst5.nino.mth.81-10.ascii",header=TRUE)
  ENOS1982_2014 <-subset(ENOS, YR>1968 & YR<2015)   # subset from 1970-2014
  ENOS1982_2014$Date<-paste(ENOS1982_2014$YR, ENOS1982_2014$MON, sep = "-" )
  ENOS1982_2014$Date<-paste(ENOS1982_2014$Date, "15", sep = "-" )
  ENOS1982_2014$Date<-as.Date(ENOS1982_2014$Date, format = "%Y-%m-%d")
  ENOS1982_2014$Anomaly <- factor(ifelse(ENOS1982_2014$ANOM.3 >=0, "Positive", "Negative"))
  
  interp <- approx(ENOS1982_2014$Date, ENOS1982_2014$ANOM.3, n=10000)
  DataENSO <- data.frame(Date=interp$x, Anomaly3=interp$y)
  DataENSO$Anomaly <- factor(ifelse(DataENSO$Anomaly3>=0, "Positive", "Negative"))
      
ENSO<-ggplot(data=DataENSO, aes(x=Date, y=Anomaly3))+
  scale_x_continuous("",expand=c(0, 0))+
  scale_y_continuous("SST anomaly (ºC) in Niño.3",expand=c(0, 0))+
  geom_area(aes(fill=Anomaly), alpha=0.7)+
  scale_fill_manual(values=ENSO_colours) + labs(title="a")+
  geom_hline(yintercept=0, linetype="solid", color="black")+
  geom_line(position = "identity", size=0.3)+
  theme(axis.text.x=element_blank(),
        axis.title.x=element_blank(),
        axis.ticks.x = element_blank())+
  annotate("text", x = 1900, y = 2.0, label = "El Niño \n72-73", size=3)+
  annotate("text", x = 5600, y = 2.0, label = "El Niño \n82-83", size=3)+
  annotate("text", x = 11100, y = 2.0, label = "El Niño \n97-98", size=3)
ENSO

# ENSOb<-ggplot(data=ENOS1982_2014, aes(x=Date, y=ANOM.3))+
#   geom_hline(yintercept=0, linetype="solid", color="black")+
#   scale_x_date("", limit=c(as.Date("1969-01-15"), as.Date("2014-12-31")),
#                breaks= "12 months",
#                      expand=c(0, 0))+
#   scale_y_continuous("SST anomaly (ºC) in Niño.3",expand=c(0, 0))+
#   geom_area(aes(fill=Anomaly))+
#   geom_line(position = "identity")+
#   scale_fill_manual(values=ENSO_colours) + labs(title="(a)")
# ENSOb

2. DHW and coral cover change

# Load packages 
library(tidyverse)
library(dplyr)
library(data.table)
library(RColorBrewer)

2.1 Yearly maximum DHW experienced by the best-represented coral reefs sites from 1982 to 2014

  • Select and transform DHW data
# Select DHW data file: maxDHW.csv
  dhw <- read.csv("Data/maxDHW.csv", header=TRUE)

# Exclude years 1981, and 2017-2019
  dhw <- dhw[!dhw$Year == 1981, ]

# Create groups for plot
  ts.Cano_island <- subset(dhw, dhw$Site == "Cano_island")     
  ts.Uva_island <- subset(dhw, dhw$Site == "Uva_Island")          
  ts.Gorgona_island <- subset(dhw, dhw$Site == "Gorgona_Island")          
  ts.Darwin_island <- subset(dhw, dhw$Site == "Darwin_North_Anchorage")          
  ts.Bahia_Banderas <- subset(dhw, dhw$Site == "Bahia_Banderas")

FIGURE 3a

  • Plot max DHW in each location
# Year, maxDHW
  par(oma = c(1,1,1,1)) 
  par(mar = c(4,4,1,1), cex.axis = 1,  las = 1)  
  par(lwd = 1)   
  
  plot(ts.Cano_island$Year,ts.Cano_island$maxDHW,   # Create a white false line with real data
       xlab="Years", ylab="Degree heating weeks (°C-weeks)",  cex=1.5, cex.lab=0.6, col="white",axes = FALSE,
       ylim = c(0, 28), lwd = 1.5) 
  par(new=T)
  rect(1982, 0, 1983.5, 28, col="#fee8c8", border ="#fee8c8")           # Strong ENSO years
  par(new=T)
  rect(1997, 0, 1998.5, 28, col="#fee8c8", border ="#fee8c8")
  par(new=T)
  plot(ts.Cano_island$Year, ts.Cano_island$maxDHW,type="l",              # Thermally stable
  xlab="", ylab="", lwd = 1.5, col="#8c2d04",axes = FALSE,
       ylim = c(0, 28))
  par(new=T)
  plot(ts.Uva_island$Year, ts.Uva_island$maxDHW, type="l",                # Thermally stable
       xlab="", ylab="", lwd = 1.5, col="#feb24c",axes = FALSE,
       ylim = c(0, 28))
  par(new=T)
  plot(ts.Gorgona_island$Year,ts.Gorgona_island$maxDHW, type="l",         # Tropical upwelling
       xlab="", ylab="",lwd = 1.5, col="#fb6a4a", axes = FALSE,
       ylim = c(0, 28))
  par(new=T)
  plot(ts.Darwin_island$Year,ts.Darwin_island$maxDHW, type="l",           # Galapagos
       xlab="", ylab="", lwd = 1.5, col="#9e9ac8",axes = FALSE,
       ylim = c(0, 28))
  par(new=T)
  plot(ts.Bahia_Banderas$Year, ts.Bahia_Banderas$maxDHW, type="l",        # Seasonal
       xlab="", ylab="",lwd = 1.5,col="#a6d96a", axes = FALSE,
       ylim = c(0, 28))
  y <- c(0,4,8,12,16,20,24,28)                                            # Axis
  axis(2,lwd = 0.5, at = y,cex.axis=0.6)
  box(lwd = 0.5)
  x <-c(1982:2014)
  YearsR = as.character(c(1982:2014))
  axis(1,  at = x, labels = YearsR, lwd = 0.4, cex.axis=0.5, las = 2) 
  abline(h=4, col="#525252", lty=2)
  abline(h=8, col="#525252", lty=2)

# Legend
  sites <- c("Thermally stable - Caño Island","Thermally stable - Uva Island",
             "Tropical upwelling - Gorgona Island","Equatorial Upwelling - Darwin Island", 
             "Seasonal - Bahía Banderas")
  par(xpd=TRUE)     
  legend(x = 2000, y = 20, 
         xpd = NA,  # or TRUE 
         inset=c(-0.5,0),
         cex = 0.4,
         legend = sites,         # box categories 
         col = c("#8c2d04","#feb24c","#fb6a4a","#9e9ac8","#a6d96a"),
         lty= 1,
         lwd = 1.5,
         x.intersp = 0.5,
         y.intersp = 2, 
         bty="n")

2.2 Relationship of coral cover rate of change and thermal stress

  • The DHW index measures the accumulated thermal stress above a bleaching threshold for 12 consecutive weeks.
  • At each site, the maximum monthly mean sea surface # temperature (MMM) is calculated from a long-term satellite climatology dataset.
  • Any temperature surpassing the MMM by 1ºC or more is added for 12 consecutive weeks,
  • because temperatures ~1ºC above the MMM are known to cause coral bleaching (Liu et al., 2006).
  • For example, if the water temperature during four consecutive weeks exceeds the MMM in 1.3°C, # 1.1°C, 1.2°C, and 1.4°C, then,
  • the site accumulates 5°C-weeks when these anomalies are summed.

We tested the hypothesis that the strength of the Log_annual_rate_of_change (response variable) is a function of the maxDHW experienced. We had measured the rate of change from four Thermal regimes, and fitted Thermal regimes as a random intercept, and estimate a common slope (change in the rate of change) for maxDHW across all Thermal regimes by fitting it as a fixed effect.

## Load data
  roc.cover <- read.csv("Data/RoC_DHW.csv", header=TRUE)  ### Select file:  RoC_DWH.csv

# Load packages
  library(lme4)
  library(lmerTest)
  library(MuMIn)

Model 6: coral cover change and DHW

model6 <- lmerTest::lmer(Log_annual_rate_of_change ~ maxDHW + (1|Thermal_regime), data=roc.cover)
step(model6)
## Backward reduced random-effect table:
## 
##                      Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)
## <none>                             4 -17.738 43.476                     
## (1 | Thermal_regime)          0    3 -27.772 61.545 20.069  1  7.469e-06
##                         
## <none>                  
## (1 | Thermal_regime) ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##        Eliminated  Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## maxDHW          0 0.65054 0.65054     1 128.59   9.873 0.002082 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## Log_annual_rate_of_change ~ maxDHW + (1 | Thermal_regime)
summary(model6) 
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Log_annual_rate_of_change ~ maxDHW + (1 | Thermal_regime)
##    Data: roc.cover
## 
## REML criterion at convergence: 35.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.7962 -0.1472  0.1167  0.4152  1.8621 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Thermal_regime (Intercept) 0.02501  0.1581  
##  Residual                   0.06589  0.2567  
## Number of obs: 133, groups:  Thermal_regime, 4
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)   
## (Intercept)  -0.085763   0.084256   3.155381  -1.018  0.38032   
## maxDHW       -0.016316   0.005193 128.594887  -3.142  0.00208 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr)
## maxDHW -0.165
r.squaredGLMM(model6)   #  returns the marginal and the conditional R2
##             R2m       R2c
## [1,] 0.05249512 0.3131505
# marginal R2  describes the proportion of variance explained by the fixed factor(s) alone:
# conditional R2 describes the proportion of variance explained by both the fixed and random factors
# Nakagawa, S. & Schielzeth, H. (2013).

FIGURE 3b: Relationship of coral cover rate of change and thermal stress

# Load packages
  library(tidyverse)
  library(dplyr)
  library(data.table)
  library(RColorBrewer)
 
# Create groups for plot
  Seasonal           = subset(roc.cover, roc.cover$Thermal_regime == "Seasonal")          
  Thermally_Stable   = subset(roc.cover, roc.cover$Thermal_regime == "Thermally_Stable")  
  Tropical_Upwelling = subset(roc.cover, roc.cover$Thermal_regime == "Upwelling") 
  Galapagos          = subset(roc.cover, roc.cover$Thermal_regime == "Galapagos_Islands") 

# Plot
  par(oma = c(1,1,1,1)) 
  par(mar = c(4,4,1,1), cex.axis = 1,  las = 1) 
  par(lwd = 1)
  plot (Seasonal$maxDHW, Seasonal$Log_annual_rate_of_change,
        pch = 21, lwd = 0.5, bg = "#a6d96a93", col="black", axes = FALSE,
        cex=1.5,cex.lab=0.6,
        xlim = c(0,30), ylim = c(-1.75, 1),
        xlab = "Degree heating weeks (°C-weeks)",
        ylab = "Annual rate of change in coral cover")
  par(new=T); 
  plot (Thermally_Stable$maxDHW, Thermally_Stable$Log_annual_rate_of_change,
        xlab="", ylab="", pch = 21, lwd = 0.5, col="black", bg = "#feb24c93", axes = FALSE,
        cex=1.5,cex.lab=0.6,
        xlim = c(0,30), ylim = c(-1.75, 1))
  par(new=T); 
  plot (Tropical_Upwelling$maxDHW, Tropical_Upwelling$Log_annual_rate_of_change,
        xlab="", ylab="", pch = 21, lwd = 0.5, bg = "#fb6a4a93",col="black", axes = FALSE,
        cex=1.5,cex.lab=0.6,
        xlim = c(0,30), ylim = c(-1.75, 1))
  par(new=T); 
  plot (Galapagos$maxDHW, Galapagos$Log_annual_rate_of_change,
        xlab="", ylab="", pch = 21, lwd = 0.5, bg = "#9e9ac893",col="black", axes = FALSE,
        cex=1.5,cex.lab=0.6,
        xlim = c(0,30), ylim = c(-1.75, 1))
  box(lwd = 0.5)
  x <-c(0,5,10,15,20,25,30)
  axis(1, lwd = 0.5, at = x, cex.axis = 0.5)
  y <- c(-1.5,-1.0,-0.5,0,0.5,1)
  axis(2,lwd = 0.5, at = y, cex.axis = 0.5)
  abline(v=4, col="#525252", lty=2); abline(v=8, col="#525252", lty=2); abline(h=0, col="#525252", lty=2); 
  # Axis text I, II, III, IV
  text(3.5, 0.3, "II", cex = 0.6); text(5, 0.3, "I", cex = 0.6)
  text(3.5, - 0.3, "III", cex = 0.6); text(5, -0.3, "IV", cex = 0.6)
  
  par(xpd=TRUE) 
  legend(x = 10, y = -0.5, 
         xpd = NA,   
         inset=c(-0.5,0),
         legend = c("Equatorial Upwelling", "Thermally stable", "Tropical upwelling","Seasonal"),
         pch=21,
         pt.bg = c("#9e9ac8","#fb6a4a","#fb6a4a","#a6d96a"),
         x.intersp = 0.5,
         y.intersp = 2, 
         bty="n")

FIGURE 3 Heat stress and its effect on coral cover across the ETP. (a) The sequence of yearly maximum DHW experienced by the best-represented individual coral reefs sites from 1982 to 2014. When cumulative thermal stress reach 4°C-weeks and 8°C-weeks or higher (horizontal dotted lines), bleaching is likely, as well as coral mortality from thermal stress (Liu et al., 2006), respectively. Although the thresholds were developed in seasonal systems, and are not universally applicable for equatorial regions, they still useful for separating responses into groups. (b) The association between the maximum DHW and the coral cover change at each of four thermal regimes (dot colors as Fig. 2). Each point (n = 133) represents the DHW and the coral cover annual rate of change per individual coral reef site. The horizontal dotted line divides the positive (> 0) and the negative (< 0) rate of change caused by coral growth or mortality, respectively. The vertical dotted line indicates the DHW thresholds of 4°C-weeks and 8°C-weeks, respectively. The vertical and horizontal dotted lines create four response quadrants: I) heat stress and positive coral cover change (resilient response), II) no heat stress and positive coral cover change, III) no heat stress and negative coral cover change, and IV) heat stress and negative coral cover change.

Descriptive statistics heat stress

No stress

no.stress <- subset(roc.cover, maxDHW < 4)                        # subset of DHW < 4
length(no.stress$maxDHW)/length(roc.cover$maxDHW)                 # ratio no stress: total              
## [1] 0.8496241
1- length(no.stress$maxDHW)/length(roc.cover$maxDHW)              # ratio stress: total              
## [1] 0.1503759
positive.rate <- subset(no.stress, Log_annual_rate_of_change > 0) # subset of Rate > 0
length(positive.rate$Log_annual_rate_of_change)                     # Number of observations in 
## [1] 48
negative.rate <- subset(no.stress, Log_annual_rate_of_change < 0)   #subset of Rate > 0
length(negative.rate$Log_annual_rate_of_change)                     # Number of observations in 
## [1] 58

No stress seasonal

no.stress.seasonal <- subset(no.stress, Thermal_regime == "Seasonal")   
positive.rate <- subset(no.stress.seasonal, Log_annual_rate_of_change > 0)  
negative.rate <- subset(no.stress.seasonal, Log_annual_rate_of_change < 0)  
length(positive.rate$Log_annual_rate_of_change)/length(negative.rate$Log_annual_rate_of_change)   
## [1] 0.07142857

No stress upwelling

no.stress.upwelling <- subset(no.stress, Thermal_regime == "Upwelling")
positive.rate <- subset(no.stress.upwelling, Log_annual_rate_of_change > 0)  
negative.rate <- subset(no.stress.upwelling, Log_annual_rate_of_change < 0)  
length(positive.rate$Log_annual_rate_of_change)/length(negative.rate$Log_annual_rate_of_change)  
## [1] 1

No stress thermally stable

no.stress.Thermally_Stable <- subset(no.stress, Thermal_regime == "Thermally_Stable")   
positive.rate <- subset(no.stress.Thermally_Stable, Log_annual_rate_of_change > 0)  
negative.rate <- subset(no.stress.Thermally_Stable, Log_annual_rate_of_change < 0)  
length(positive.rate$Log_annual_rate_of_change)/length(negative.rate$Log_annual_rate_of_change)  
## [1] 1.2

No stress Galapagos

no.stress.Galapagos_Islands <- subset(no.stress, Thermal_regime == "Galapagos_Islands") 
positive.rate <- subset(no.stress.Galapagos_Islands, Log_annual_rate_of_change > 0)  
negative.rate <- subset(no.stress.Galapagos_Islands, Log_annual_rate_of_change < 0)  
length(positive.rate$Log_annual_rate_of_change)/length(negative.rate$Log_annual_rate_of_change)  
## [1] 1

3. Supplementary figures

library(pals)
library(dplyr)
library(reshape2)
library(RColorBrewer)

FIGURE S1: Number of regionsand countries surveyed per year

# Barplot number of surveys per year
  obs_num = as.data.frame(table(cover$Year,cover$Region))    
  colnames(obs_num) = c("Year", "Region", "Freq")
  wide_obs_num = dcast(obs_num, Region~Year, value = "Freq")  
  df = data.matrix(wide_obs_num[,2:45])   
  as.table(df)     
##   1970 1971 1972 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984
## A    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## B    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## C    0    0    1    0    0    0    0    0    1    1    0    0    1    0
## D    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## E    0    0    0    0    0    0    0    0    0    1    0    1    0    4
## F    0    0    0    2    3    2    0    0    0    0    0    1    0    1
## G    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## H    3   10    2   13    4    4    3    5    2    5    2    3    5    3
## I    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## J    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## K    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## L    0    0    0    0    0    0    0    2    0    0    0    0    0    0
## M    0    0    0    0    0    0    0    0    0    0    0    1    0    0
## N    2    6    0    2    1    0    0    0    0    0    0    1    4    2
## O    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## P    0    0    0    0    0    0    0    0    0    0    0    0    0    0
##   1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998
## A    0    0    0    0    0    0    0    0    0    1    0    0    0    0
## B    0    0    5    0    0    0    0    1    0    0    0    0    0    0
## C    0    0    1    2    4    0    0    1    0    0    2    2    0    4
## D    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## E   12    1    2    0    1    0    0   11    0    2    0    6    0    1
## F    0    1    1    0    1    0    0    0    2    0    0    0    0    0
## G    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## H    4    5    2    3    1    1    0    4    1    4    4    0    1    1
## I    0    0    0    0    0    0    2    0    0    0    0    0   16   10
## J    0    0    1    1    1    0    3    1    1    0    0    0    2    2
## K    0    0    0    0    0    0    0    0    0    0    0    0    2    1
## L    0    0    0    0    0    0    0    0    0    1    1    4    4    2
## M    0    0    0    0    0    0    0    1    0    0    0    0    0    0
## N    0    1    5    1    2    0    1    1    1    1    1    1    1    1
## O    0    0    0    0    0    0    0    0    1    1    0    0    0    0
## P    0    0    0    0    0    0    0    0    0    0    0    0    0    0
##   1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012
## A    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## B    0    0    0    1    0    0    0    1    0    0    4    0    0    0
## C    7    4    3    4    6    4    3    6    3    2    3    2    2    3
## D    0    0    0    0    0    0    0    0    0    0    0    4    0    1
## E    1    0    2    0    4    2    2    0    1    0    2    0    0    0
## F    0    0    1    1    0    1    3    6    2    2    0    0    1    0
## G    0    0    0    0    0    0    0    0    0    0    1    0    0    0
## H    0    1    4   14   12    1    1    1    1    0    1    2    0    0
## I    0    0    0    6    3    0    1    1    4    0   17    6    2    0
## J    0    0    0    1    2    0    0    0    0    5    7    0    0    1
## K    1    5    5    3    2    0    0    0    0    0    2    0    0    0
## L    0    0    2    0    1    0    1    3    2    0   17    3    1    1
## M    0    0    0    0    0    1    2   14    8    0    0    0    0    1
## N    1    1    1    1    1    1    0    0    0    0    1    0    0    0
## O    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## P    0    0    0    0    0    2    0    0    0    0    0    0    0    0
##   2013 2014
## A    0    0
## B    1    1
## C    2    1
## D    4    2
## E    0    0
## F    0    0
## G    0    0
## H    0    0
## I    2    2
## J    0    0
## K    0    0
## L    1    3
## M    0    0
## N    0    0
## O    0    0
## P    0    0
  df[ df>1 ] <- 1 
  levels(cover$Region) # How many regions? = how many colours
##  [1] "Clipperton"                 "Cocos_Islands"             
##  [3] "Colombia"                   "Costa_Rica_Central"        
##  [5] "Costa_Rica_Southern"        "Eastern_Galapagos_Islands" 
##  [7] "El_Salvador"                "Gulf_of_Chiriqui"          
##  [9] "Mexican_Central"            "Mexican_Northern"          
## [11] "Mexican_Southern"           "Nicaragua_Papagayo_Zone"   
## [13] "Northern_Galapagos_Islands" "Panama_Bight"              
## [15] "Revillagigedo_Islands"      "Western_Galapagos_Islands"
  Regions <- c(levels(obs_num$Region))
  color.regions = warmcool(17)
  
  # Figure 2a
  par(mfcol=c(2,1), cex.main = 1.5, mar = c(4,4,1,1), cex.lab = 1.5, cex.axis = 1,  las = 1)
  barplot(df,
          ylim=c(0,15),
          col = color.regions,   
          ylab="Number of regions surveyed per year",  
          xlab="Year",
          #lwd = 1,          
          las = 2,           
          axis.lty = 1, 
          cex.names = 0.4,   
          cex.axis = 0.8)
# Legend as a box
  legend( #"bottomright",    
     x = 2, y = 16, 
     xpd = NA,  # or TRUE 
     legend = Regions,          
     fill = color.regions[1:17],   
     inset=c(-0.5,0),
     cex = 0.35,
     text.font=0.5,            
     x.intersp = -0.5,
     y.intersp = 1,
     lwd = 0.1,
     bty="n")

# Figure 2b
obs_num = as.data.frame(table(cover$Year,cover$Country))    # Table per country
colnames(obs_num) = c("Year", "Country", "Freq")
wide_obs_num = dcast(obs_num, Country~Year, value = "Freq")   
dfc = data.matrix(wide_obs_num[,2:45])   
as.table(dfc)       
##   1970 1971 1972 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984
## A    0    0    1    0    0    0    0    0    1    1    0    0    1    0
## B    0    0    2    9    0    0    0    2    0    2    0    1    0    4
## C    0    0    0    2    3    2    0    0    0    0    0    2    0    1
## D    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## E    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## F    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## G    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## H    5   16    0    6    5    4    3    5    2    4    2    4    9    5
##   1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998
## A    0    0    1    2    4    0    0    1    0    0    2    2    0    4
## B   12    1    7    0    1    0    0   15    0    6    4   10    4    3
## C    0    1    1    0    1    0    0    1    2    0    0    0    0    0
## D    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## E    0    0    0    0    0    0    0    0    0    1    0    0    0    0
## F    0    0    1    1    1    0    5    1    2    1    0    0   20   13
## G    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## H    4    6    7    4    3    1    1    2    2    2    2    1    2    2
##   1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012
## A    7    4    3    4    6    4    3    6    3    2    3    2    2    3
## B    1    0    7    1    5    2    3    4    3    0   12    7    1    2
## C    0    0    1    1    0    4    5   20   10    2    0    0    1    1
## D    0    0    0    0    0    0    0    0    0    0    1    0    0    0
## E    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## F    1    5    5   10    7    0    1    1    4    5   26    6    2    1
## G    0    0    0    0    0    0    0    0    0    0   11    0    0    0
## H    1    2    2   15   13    2    1    1    1    0    2    2    0    0
##   2013 2014
## A    2    1
## B    6    6
## C    0    0
## D    0    0
## E    0    0
## F    2    2
## G    0    0
## H    0    0
dfc[ dfc>1 ] <- 1 # replace all values greater than one with one
color.countries <- brewer.pal(9, "Oranges")
countries = c(levels(obs_num$Country))
barplot(dfc,
        col = color.countries,   
        ylim=c(0,10),
        ylab="Number of countries surveyed per year",  
        xlab="Year",
        las = 2,           
        axis.lty = 1,      
        cex.names = 0.4,   
        cex.axis = 0.8)
# Legend as a box
legend( #"bottomright",    
   x = 2, y = 10, 
   xpd = NA,  # or TRUE 
   legend = countries,         
   fill = color.countries[1:9],   
   inset=c(-0.5,0),
   cex = 0.35,
   text.font=1,            
   x.intersp = -0.5,
   y.intersp = 1,
   lwd = 0.1,
   bty="n")

FIGURE S1 Regions and countries surveyed per year. (a) Number of regions surveyed per year in the ETP. Since 1970, the number of surveys in all regions increased by year up to a maximum of 10 regions surveyed in 2009. (b) Number of countries surveyed per year (Colombia, Costa Rica, Ecuador, El Salvador, France, Mexico, Nicaragua and Panama). Since 1970, the number of surveys in all the countries reached a maximum of 6 countries in 2009.

FIGURE S2: Variation in the live coral cover for 1970–2014 per thermal regime

Temporal variation in the live coral cover for 1970–2014.

    1. Tropical upwelling regime.
    1. Seasonal regime.
    1. Thermally stable regime.
    1. Equatorial upwelling regime.

The smooth line represents the best fit to the time series with a 95% confidence level interval and a span = 0.3.

library(ggplot2)
library(lemon)

aggr.region_S2<-aggr.region
aggr.region_S2$Thermal_regime<-factor(aggr.region_S2$Thermal_regime, 
                                  levels=c("Tropical upwelling", "Seasonal", 
                                           "Thermally stable", "Equatorial upwelling"),
                                  labels = c ("a", "b", "c", "d"))
                                   

FigureS2<-ggplot(aggr.region_S2, aes(x=Year, y=Coral_cover)) + 
  theme_bw() + theme(panel.border = element_blank(),   # set ggplot to white background
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(), 
                              axis.line = element_line(),
                              strip.background = element_blank(),
                              strip.text = element_text(face="bold",hjust = 0, vjust = -1))+
  geom_point()+ geom_smooth (span = 0.3) +
  scale_x_continuous(limits=c(1970, 2014), 
                     breaks = seq(1970, 2014, by=5),
                     expand = c(0.01, 0.01), "Year") +
  scale_y_continuous(limits=c(-19, 91), expand = c(0, 0), "Live coral cover (%)") +
  facet_wrap(Thermal_regime~., ncol=1, scales = "free_x", strip.position ="top")

FigureS2

#ggsave(file="FigureS2.svg", plot=FigureS2, width=6.69, height=6.69)

FIGURE S2 Temporal variation in the live coral cover for 1970–2014. (a) Tropical upwelling regime (n = 75). (b) Seasonal regime (n = 28). (c) Thermally stable regime (n = 74). (d) Equatorial upwelling regime (n = 25). The smooth line represents the best fit to the time series with a 95% confidence level interval and a span = 0.3.

FigureS2b<-ggplot(cover, aes(x=Year, y=Coral_cover)) + 
  theme_bw() + theme(panel.border = element_blank(),   # set ggplot to white background
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(), 
                              axis.line = element_line(),
                              strip.background = element_blank(),
                              strip.text = element_text(face="bold",hjust = 0, vjust = -1))+
  geom_point(aes(colour=Region)) +
  geom_smooth (span = 0.3, colour="gray") +
  scale_x_continuous(limits=c(1970, 2014), 
                     breaks = seq(1970, 2014, by=5),
                     expand = c(0.01, 0.01), "Year") +
  scale_y_continuous(limits=c(-19, 91), expand = c(0, 0), "Live coral cover (%)") +
  facet_wrap(Thermal_regime~., ncol=1, scales = "free_x", strip.position ="top") +
  theme(legend.position = "bottom", legend.direction = "horizontal" )

FigureS2b

4. Revision figures

Data only for four best monitoring sites

We found in published studies four comprehensive monitoring datasets for Uva Island, Gorgona Island, Costa Rica, and the Eastern and Northern Galapagos Islands, which correspond to well recognized long-term monitoring programs

# Create groups
Uva_Island <- subset(aggr.site, Site == "Uva_Island")
Gorgona <-subset(aggr.site, Site == "La_Azufrada")
Cano_island <-subset(aggr.site, Site == "Cano_island")
EGI <-subset(aggr.site, Region == "Eastern_Galapagos_Islands")
NGI <-subset(aggr.site, Region == "Northern_Galapagos_Islands")

Monitoring_sites <- rbind(Uva_Island,Gorgona,Cano_island,EGI,NGI)

Model 7: Only monitored sites

model7 <- lme(Coral_cover ~ -1 + Thermal_regime * Year, random = ~1|Region, data=Monitoring_sites) 
summary(model7)
## Linear mixed-effects model fit by REML
##  Data: Monitoring_sites 
##        AIC      BIC    logLik
##   978.1405 999.7444 -481.0703
## 
## Random effects:
##  Formula: ~1 | Region
##         (Intercept) Residual
## StdDev:    14.54312 15.90812
## 
## Fixed effects: Coral_cover ~ -1 + Thermal_regime * Year 
##                                             Value Std.Error  DF    t-value
## Thermal_regimeThermally stable          -565.5374  449.5031   2 -1.2581391
## Thermal_regimeTropical upwelling         262.8581  757.1479   2  0.3471688
## Thermal_regimeEquatorial upwelling       373.7671  418.4179   2  0.8932865
## Year                                       0.2928    0.2257 109  1.2971984
## Thermal_regimeTropical upwelling:Year     -0.3987    0.4404 109 -0.9054869
## Thermal_regimeEquatorial upwelling:Year   -0.4645    0.3078 109 -1.5092165
##                                         p-value
## Thermal_regimeThermally stable           0.3353
## Thermal_regimeTropical upwelling         0.7616
## Thermal_regimeEquatorial upwelling       0.4660
## Year                                     0.1973
## Thermal_regimeTropical upwelling:Year    0.3672
## Thermal_regimeEquatorial upwelling:Year  0.1341
##  Correlation: 
##                                         Thr_Ts Thr_Tu Thr_Eu Year   T_Tu:Y
## Thermal_regimeTropical upwelling         0.000                            
## Thermal_regimeEquatorial upwelling       0.000  0.000                     
## Year                                    -1.000  0.000  0.000              
## Thermal_regimeTropical upwelling:Year    0.512 -0.858  0.000 -0.513       
## Thermal_regimeEquatorial upwelling:Year  0.733  0.000 -0.680 -0.733  0.376
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.73438942 -0.56940133 -0.09515832  0.50988895  3.13690545 
## 
## Number of Observations: 116
## Number of Groups: 5
anova(model7)
##                     numDF denDF  F-value p-value
## Thermal_regime          3     2 7.535454  0.1194
## Year                    1   109 0.023766  0.8778
## Thermal_regime:Year     2   109 1.205466  0.3035
plot(ranef(model7))    

plot(model7)           

resnorm7 <- resid(model7)
hist(resnorm7, xlab = "Residuals", main = "") 

coef.m7 <- as.data.frame(coef(summary(model7))) 
plot(model7, resid(., type = "p") ~ fitted(.) | Region, abline = 0)

plot(model7, Coral_cover ~ fitted(.) | Region, abline = c(0,1))

plot(model7, Coral_cover ~ fitted(.) | Thermal_regime, abline = c(0,1))

FIGURE S3: Data only for four best monitoring sites

FigureS3 <-ggplot(Monitoring_sites, aes(x=Year, y=Coral_cover)) +
  geom_boxplot(aes(group=Year), outlier.shape = NA) +
  stat_boxplot(aes(group=Year), geom = 'errorbar')+
  geom_smooth(method = lm, se=FALSE, linetype = "dashed", colour="red")+
  geom_smooth(span = 0.2, se=T, colour="darkgray")+ 
  scale_y_continuous("Live coral cover (%)", limits = c(-15, 95), 
                     breaks = seq(0, 90, by=20), expand = c(0,0))+
  scale_x_continuous("", limits = c(1969, 2015),
                     breaks = seq(1970, 2014, by=2), expand = c(0,0))+
  annotate("rect", xmin = 1982, xmax = 1983, ymin = 0, ymax = 80,  alpha = .2, fill="red")+
  annotate("rect", xmin = 1997, xmax = 1998, ymin = 0, ymax = 80,  alpha = .2, fill="red")+
  annotate("rect", xmin = 1972, xmax = 1973, ymin = 0, ymax = 80,  alpha = .2, fill="red")

FigureS3

FigureS3b<- FigureS3 + geom_point(aes(colour=Location), alpha=0.4) +
  theme(legend.position = "bottom")
FigureS3b

#ggsave(file="FigS3.png", plot=FigureS3, width=6.69, height=4)
#ggsave(file="FigS3b.png", plot=FigureS3b, width=6.69, height=4.69)

Data excluding one thermal regime from model 1

Model 1: All years (1970_2014) removing each termal regime

  • Tropical upwelling
Model1_NoEquatorial_Upwelling <- ggplot(data=subset(aggr.region, Thermal_regime!="Tropical upwelling"), aes(x=Year, y=Coral_cover)) +
  geom_boxplot(aes(group=Year), outlier.shape = NA) +
  stat_boxplot(aes(group=Year), geom = 'errorbar')+
  geom_smooth(method = lm, se=FALSE, linetype = "dashed", colour="red")+
  geom_smooth(span = 0.2, se=T, colour="darkgray")+ 
  #geom_point(aes(colour=Thermal_regime), alpha=0.5) +
  scale_y_continuous("Live coral cover (%)", limits = c(-5, 90), 
                     breaks = seq(0, 90, by=20), expand = c(0,0))+
  scale_x_continuous("", limits = c(1969, 2015),
                     breaks = seq(1970, 2014, by=2), expand = c(0,0))+
  scale_color_manual(values=my_colours) + labs(title="a")+
  #annotate("text", x = 1983, y = 70, label = "*", size=8)+
  annotate("rect", xmin = 1982, xmax = 1983, ymin = 0, ymax = 80,  alpha = .2, fill="red")+
  annotate("rect", xmin = 1997, xmax = 1998, ymin = 0, ymax = 80,  alpha = .2, fill="red")+
  annotate("rect", xmin = 1972, xmax = 1973, ymin = 0, ymax = 80,  alpha = .2, fill="red")
#Model1_NoEquatorial_Upwelling
  • Seasonal
Model1_NoSeasonal <- ggplot(data=subset(aggr.region, Thermal_regime!="Seasonal"), aes(x=Year, y=Coral_cover)) +
  geom_boxplot(aes(group=Year), outlier.shape = NA) +
  stat_boxplot(aes(group=Year), geom = 'errorbar')+
  geom_smooth(method = lm, se=FALSE, linetype = "dashed", colour="red")+
  geom_smooth(span = 0.2, se=T, colour="darkgray")+ 
  #geom_point(aes(colour=Thermal_regime), alpha=0.5) +
  scale_y_continuous("Live coral cover (%)", limits = c(-5, 90),
                     breaks = seq(0, 90, by=20), expand = c(0,0))+
  scale_x_continuous("", limits = c(1969, 2015),
                     breaks = seq(1970, 2014, by=2), expand = c(0,0))+
  scale_color_manual(values=my_colours) + labs(title="b")+
  #annotate("text", x = 1983, y = 70, label = "*", size=8)+
  annotate("rect", xmin = 1982, xmax = 1983, ymin = 0, ymax = 80,  alpha = .2, fill="red")+
  annotate("rect", xmin = 1997, xmax = 1998, ymin = 0, ymax = 80,  alpha = .2, fill="red")+
  annotate("rect", xmin = 1972, xmax = 1973, ymin = 0, ymax = 80,  alpha = .2, fill="red")
#Model1_NoSeasonal
  • Thermally stable
Model1_NoThermallyStable <- ggplot(data=subset(aggr.region, Thermal_regime!="Thermally stable"), aes(x=Year, y=Coral_cover)) +
  geom_boxplot(aes(group=Year), outlier.shape = NA) +
  stat_boxplot(aes(group=Year), geom = 'errorbar')+
  geom_smooth(method = lm, se=FALSE, linetype = "dashed", colour="red")+
  geom_smooth(span = 0.2, se=T, colour="darkgray")+ 
  #geom_point(aes(colour=Thermal_regime), alpha=0.5) +
  scale_y_continuous("Live coral cover (%)", limits = c(-5, 90),
                     breaks = seq(0, 90, by=20), expand = c(0,0))+
  scale_x_continuous("", limits = c(1969, 2015),
                     breaks = seq(1970, 2014, by=2), expand = c(0,0))+
  scale_color_manual(values=my_colours) + labs(title="c")+
  #annotate("text", x = 1983, y = 70, label = "*", size=8)+
  annotate("rect", xmin = 1982, xmax = 1983, ymin = 0, ymax = 80,  alpha = .2, fill="red")+
  annotate("rect", xmin = 1997, xmax = 1998, ymin = 0, ymax = 80,  alpha = .2, fill="red")+
  annotate("rect", xmin = 1972, xmax = 1973, ymin = 0, ymax = 80,  alpha = .2, fill="red")
# Model1_NoThermallyStable
  • Equatorial upwelling
Model1_No_EqUpwelling <- ggplot(data=subset(aggr.region, Thermal_regime!="Equatorial upwelling"), aes(x=Year, y=Coral_cover)) +
  geom_boxplot(aes(group=Year), outlier.shape = NA) +
  stat_boxplot(aes(group=Year), geom = 'errorbar')+
  geom_smooth(method = lm, se=FALSE, linetype = "dashed", colour="red")+
  geom_smooth(span = 0.2, se=T, colour="darkgray")+ 
  #geom_point(aes(colour=Thermal_regime), alpha=0.5) +
  scale_y_continuous("Live coral cover (%)", limits = c(-5, 90), 
                     breaks = seq(0, 90, by=20), expand = c(0,0))+
  scale_x_continuous("", limits = c(1969, 2015),
                     breaks = seq(1970, 2014, by=2), expand = c(0,0))+
  scale_color_manual(values=my_colours) + labs(title="d")+
  #annotate("text", x = 1983, y = 70, label = "*", size=8)+
  annotate("rect", xmin = 1982, xmax = 1983, ymin = 0, ymax = 80,  alpha = .2, fill="red")+
  annotate("rect", xmin = 1997, xmax = 1998, ymin = 0, ymax = 80,  alpha = .2, fill="red")+
  annotate("rect", xmin = 1972, xmax = 1973, ymin = 0, ymax = 80,  alpha = .2, fill="red")
# Model1_No_EqUpwelling

FIGURE S4: Effect of thermal regime removal/addition

FigureS4<-grid.arrange(Model1_NoEquatorial_Upwelling, Model1_NoSeasonal, 
                       Model1_NoThermallyStable, Model1_No_EqUpwelling,
                       ncol = 2)

#ggsave(file="FigS5.png", plot=FigureS5, width=6.69, height=6.69)
#ggsave(file="FigS5.pdf", plot=FigureS5, width=6.69, height=6.69)

Figure S4: Coral cover trends excluding (a) Thermally stable, (b) Tropical Upwelling, (c) Equatorial upwelling, and (d) Seasonal thermal regimes

Add CI based on boostraping to fig 2?

To use “bootMer” models need to be run with lme4 instead nlme. Please check accuracy of the new (b) models

Model 1b: All years 1970_2014 with lme4

  # 1. Build lme4 model:
   
    model1b <- lmer(
      Coral_cover ~ Thermal_regime * Year + 
        (1|Region), data=aggr.region, REML = FALSE)
  
    summary(model1b)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: Coral_cover ~ Thermal_regime * Year + (1 | Region)
##    Data: aggr.region
## 
##      AIC      BIC   logLik deviance df.resid 
##   1746.3   1779.4   -863.2   1726.3      192 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.32243 -0.72526  0.00017  0.57206  3.08905 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Region   (Intercept)  96.1     9.803  
##  Residual             268.7    16.393  
## Number of obs: 202, groups:  Region, 16
## 
## Fixed effects:
##                                           Estimate Std. Error         df
## (Intercept)                             -425.34274  371.24038  200.81720
## Thermal_regimeTropical upwelling         155.93891  520.89698  201.08927
## Thermal_regimeEquatorial upwelling       661.72680  677.59229  195.76846
## Thermal_regimeSeasonal                  1436.94164  901.75966  201.77989
## Year                                       0.22505    0.18575  200.56729
## Thermal_regimeTropical upwelling:Year     -0.07308    0.26052  200.89738
## Thermal_regimeEquatorial upwelling:Year   -0.33411    0.33909  195.36095
## Thermal_regimeSeasonal:Year               -0.71540    0.45058  201.68330
##                                         t value Pr(>|t|)
## (Intercept)                              -1.146    0.253
## Thermal_regimeTropical upwelling          0.299    0.765
## Thermal_regimeEquatorial upwelling        0.977    0.330
## Thermal_regimeSeasonal                    1.593    0.113
## Year                                      1.212    0.227
## Thermal_regimeTropical upwelling:Year    -0.280    0.779
## Thermal_regimeEquatorial upwelling:Year  -0.985    0.326
## Thermal_regimeSeasonal:Year              -1.588    0.114
## 
## Correlation of Fixed Effects:
##             (Intr) Thr_Tu Thr_Eu Thrm_S Year   T_Tu:Y T_Eu:Y
## Thrml_rgmTu -0.713                                          
## Thrml_rgmEu -0.548  0.390                                   
## Thrml_rgmSs -0.397  0.283  0.218                            
## Year        -1.000  0.713  0.548  0.397                     
## Thrml_rTu:Y  0.713 -1.000 -0.391 -0.283 -0.713              
## Thrml_rEu:Y  0.548 -0.390 -1.000 -0.217 -0.548  0.391       
## Thrml_rgS:Y  0.398 -0.283 -0.218 -1.000 -0.397  0.283  0.217
    anova(model1b) # Thermal regime is not significant anymore?? 
## Type III Analysis of Variance Table with Satterthwaite's method
##                     Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## Thermal_regime      839.67 279.889     3 199.65  1.0415 0.3753
## Year                 41.27  41.273     1 199.94  0.1536 0.6956
## Thermal_regime:Year 843.71 281.237     3 199.47  1.0465 0.3731
    ranova(model1b) # Region seems the only significant effect
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## Coral_cover ~ Thermal_regime + Year + (1 | Region) + Thermal_regime:Year
##              npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>         10 -863.15 1746.3                         
## (1 | Region)    9 -870.85 1759.7 15.399  1  8.705e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    drop1(model1b)
## Single term deletions using Satterthwaite's method:
## 
## Model:
## Coral_cover ~ Thermal_regime * Year + (1 | Region)
##                     Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## Thermal_regime:Year 843.71  281.24     3 199.47  1.0465 0.3731

Not sure if this if right since the model itself is not significant 

  # 2. Predict values:
    pred1 <- predict(model1b,re.form = NA)

  #3. Bootstrap CI:
    boot1 <- bootMer(model1b, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
    std.err <- apply(boot1$t, 2, sd)
    CI.lo_1 <- pred1 - std.err*1.96
    CI.hi_1 <- pred1 + std.err*1.96

  #Plot
  Model1b_plot<- ggplot(
    aggr.region, aes(x = Year, y = Coral_cover, colour = Thermal_regime)) +
    geom_line(aes(y = pred1),size=2) +
    geom_point(aes(fill=factor(Thermal_regime)),
             shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.5) +
    geom_ribbon(aes(ymin = CI.lo_1, ymax = CI.hi_1),
                size=2, alpha = 0.03, linetype = 0) +
    scale_color_manual(values=my_colours) +
    scale_fill_manual(values=my_colours) +
    scale_y_continuous("Live coral cover (%)", 
                       limits = c(-20, 90), 
                       breaks = seq(0, 90, by=20), expand = c(0,0))+
    scale_x_continuous("", limits = c(1969, 2015),
                     breaks = seq(1970, 2014, by=2), expand = c(0,0))
  
  Model1b_plot

Model 2b: Time-interval 1 - ENSO.73_82 with lme4

Not enough data

Model 3b: Time-interval 2 - ENSO.83_97 with lme4

  # 1. Build lme4 model:
   
    model3b <- lmer(
      Coral_cover ~ Thermal_regime * Year + 
        (1|Region), data=ENSO.83_97, REML = FALSE)
  
    summary(model3b)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: Coral_cover ~ Thermal_regime * Year + (1 | Region)
##    Data: ENSO.83_97
## 
##      AIC      BIC   logLik deviance df.resid 
##    499.0    520.8   -239.5    479.0       55 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0243 -0.5395 -0.0293  0.5949  3.2102 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Region   (Intercept)  7.645   2.765   
##  Residual             86.995   9.327   
## Number of obs: 65, groups:  Region, 10
## 
## Fixed effects:
##                                           Estimate Std. Error         df
## (Intercept)                             -2558.7536   885.2036    59.4204
## Thermal_regimeTropical upwelling        -4663.3240  1215.3575    62.0348
## Thermal_regimeEquatorial upwelling       2670.5406  2851.0741    56.7512
## Thermal_regimeSeasonal                   8827.9001  2025.6205    62.6391
## Year                                        1.2929     0.4449    59.3736
## Thermal_regimeTropical upwelling:Year       2.3555     0.6105    61.9948
## Thermal_regimeEquatorial upwelling:Year    -1.3488     1.4341    56.7453
## Thermal_regimeSeasonal:Year                -4.4231     1.0170    62.5994
##                                         t value Pr(>|t|)    
## (Intercept)                              -2.891 0.005362 ** 
## Thermal_regimeTropical upwelling         -3.837 0.000294 ***
## Thermal_regimeEquatorial upwelling        0.937 0.352893    
## Thermal_regimeSeasonal                    4.358 4.97e-05 ***
## Year                                      2.906 0.005132 ** 
## Thermal_regimeTropical upwelling:Year     3.858 0.000275 ***
## Thermal_regimeEquatorial upwelling:Year  -0.941 0.350946    
## Thermal_regimeSeasonal:Year              -4.349 5.13e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Thr_Tu Thr_Eu Thrm_S Year   T_Tu:Y T_Eu:Y
## Thrml_rgmTu -0.728                                          
## Thrml_rgmEu -0.310  0.226                                   
## Thrml_rgmSs -0.430  0.313  0.133                            
## Year        -1.000  0.728  0.310  0.430                     
## Thrml_rTu:Y  0.729 -1.000 -0.226 -0.313 -0.729              
## Thrml_rEu:Y  0.310 -0.226 -1.000 -0.133 -0.310  0.226       
## Thrml_rgS:Y  0.430 -0.313 -0.133 -1.000 -0.430  0.313  0.133
    anova(model3b) # Thermal regime and its interaciton with year are significant 
## Type III Analysis of Variance Table with Satterthwaite's method
##                     Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Thermal_regime      4445.7 1481.89     3 60.690 17.0343 3.858e-08 ***
## Year                  87.0   86.98     1 60.112  0.9998    0.3214    
## Thermal_regime:Year 4462.4 1487.47     3 60.665 17.0984 3.672e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ranova(model3b) # ramdom effects are not 
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## Coral_cover ~ Thermal_regime + Year + (1 | Region) + Thermal_regime:Year
##              npar  logLik    AIC    LRT Df Pr(>Chisq)
## <none>         10 -239.52 499.03                     
## (1 | Region)    9 -240.10 498.20 1.1616  1     0.2811
    drop1(model3b)
## Single term deletions using Satterthwaite's method:
## 
## Model:
## Coral_cover ~ Thermal_regime * Year + (1 | Region)
##                     Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Thermal_regime:Year 4462.4  1487.5     3 60.665  17.098 3.672e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  # 2. Predict values:
    pred3 <- predict(model3b,re.form = NA)

  #3. Bootstrap CI:
    boot3 <- bootMer(model3b, predict, nsim = 1000, re.form = NA) # No ramndom effects
    std.err <- apply(boot3$t, 2, sd)
    CI.lo_3 <- pred3 - std.err*1.96
    CI.hi_3 <- pred3 + std.err*1.96

  #Plot
    
   Model3b_plot<- ggplot(
     ENSO.83_97, aes(x = Year, y = Coral_cover, colour = Thermal_regime))+
    geom_line(aes(y = pred3),size=2) +
    geom_point(aes(fill=factor(Thermal_regime)),
          shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.5) + 
    geom_ribbon(aes(ymin = CI.lo_3, ymax = CI.hi_3),
                size=2, alpha = 0.03, linetype = 0) +
    scale_fill_manual(values=my_colours)+
    scale_colour_manual(values=my_colours) +
    scale_x_continuous("", limits = c(1982.7, 1998.2), 
                        breaks = seq(1983, 1997, by=2), 
                        expand = c(0.02,0.02))+
    scale_y_continuous("Live coral cover (%)", 
                       limits = c(-20, 90), 
                       breaks = seq(0, 90, by=20), expand = c(0,0))
    
   Model3b_plot

  #ggsave(file="Model3b_plot.png", plot=Model3b_plot, width=4, height=4)
  #ggsave(file="Model3b_plot.pdf", plot=Model3b_plot, width=4, height=4)

Model 4b: Time-interval 3 - ENSO.98_14 with lme4

  # 1. Build lme4 model:
   
    model4b <- lmer(
      Coral_cover ~ Thermal_regime * Year + 
        (1|Region), data=ENSO.98_14, REML = FALSE)
  
    summary(model4b)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: Coral_cover ~ Thermal_regime * Year + (1 | Region)
##    Data: ENSO.98_14
## 
##      AIC      BIC   logLik deviance df.resid 
##    894.7    921.2   -437.4    874.7       94 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.97345 -0.54963 -0.01987  0.30764  2.79401 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Region   (Intercept)  81.52    9.029  
##  Residual             224.67   14.989  
## Number of obs: 104, groups:  Region, 12
## 
## Fixed effects:
##                                           Estimate Std. Error         df
## (Intercept)                             -1734.4099  1338.8526   103.8926
## Thermal_regimeTropical upwelling         4885.3253  1681.7337   103.9991
## Thermal_regimeEquatorial upwelling       -332.4797  3113.7676    96.2513
## Thermal_regimeSeasonal                   1970.4798  2034.4368   100.1488
## Year                                        0.8782     0.6677   103.9042
## Thermal_regimeTropical upwelling:Year      -2.4312     0.8385   103.9997
## Thermal_regimeEquatorial upwelling:Year     0.1641     1.5523    96.2303
## Thermal_regimeSeasonal:Year                -0.9798     1.0140   100.0511
##                                         t value Pr(>|t|)   
## (Intercept)                              -1.295  0.19804   
## Thermal_regimeTropical upwelling          2.905  0.00449 **
## Thermal_regimeEquatorial upwelling       -0.107  0.91519   
## Thermal_regimeSeasonal                    0.969  0.33510   
## Year                                      1.315  0.19132   
## Thermal_regimeTropical upwelling:Year    -2.899  0.00456 **
## Thermal_regimeEquatorial upwelling:Year   0.106  0.91605   
## Thermal_regimeSeasonal:Year              -0.966  0.33623   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Thr_Tu Thr_Eu Thrm_S Year   T_Tu:Y T_Eu:Y
## Thrml_rgmTu -0.796                                          
## Thrml_rgmEu -0.430  0.342                                   
## Thrml_rgmSs -0.652  0.519  0.280                            
## Year        -1.000  0.796  0.430  0.652                     
## Thrml_rTu:Y  0.796 -1.000 -0.342 -0.519 -0.796              
## Thrml_rEu:Y  0.430 -0.342 -1.000 -0.280 -0.430  0.343       
## Thrml_rgS:Y  0.652 -0.519 -0.280 -1.000 -0.652  0.519  0.280
    anova(model4b) # Thermal regime and its interaciton with year are significant 
## Type III Analysis of Variance Table with Satterthwaite's method
##                      Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## Thermal_regime      2257.18  752.39     3 99.562  3.3489 0.02210 *
## Year                   4.86    4.86     1 96.337  0.0217 0.88333  
## Thermal_regime:Year 2248.88  749.63     3 99.519  3.3366 0.02244 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ranova(model4b) # ramdom effects are significant
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## Coral_cover ~ Thermal_regime + Year + (1 | Region) + Thermal_regime:Year
##              npar  logLik    AIC   LRT Df Pr(>Chisq)    
## <none>         10 -437.36 894.73                        
## (1 | Region)    9 -443.08 904.15 11.42  1  0.0007265 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    drop1(model4b)
## Single term deletions using Satterthwaite's method:
## 
## Model:
## Coral_cover ~ Thermal_regime * Year + (1 | Region)
##                     Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## Thermal_regime:Year 2248.9  749.63     3 99.519  3.3366 0.02244 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  # 2. Predict values:
    pred4 <- predict(model4b,re.form = NA)

  #3. Bootstrap CI:
    boot4 <- bootMer(model4b, predict, nsim = 1000, re.form = NULL) # Include ramndom effects
    std.err <- apply(boot4$t, 2, sd)
    CI.lo_4 <- pred4 - std.err*1.96
    CI.hi_4 <- pred4 + std.err*1.96

  #Plot
   Model4b_plot<- ggplot(ENSO.98_14, aes(
     x = Year, y = Coral_cover, colour = Thermal_regime)) +
    geom_line(aes(y = pred4),size=2) +
    geom_point(aes(fill=factor(Thermal_regime)),
             shape = 21, colour = "black", size = 2, 
             stroke = 0.3, alpha=0.5) +
    geom_ribbon(aes(ymin = CI.lo_4, ymax = CI.hi_4),size=2, 
                alpha = 0.03, linetype = 0) +
    scale_y_continuous("Live coral cover (%)",
                       limits = c(-20, 90), 
                       breaks = seq(0, 90, by=20), expand = c(0,0))+
    scale_x_continuous("", limits = c(1997.7, 2014.3),
                     breaks = seq(1998, 2014, by=2), expand = c(0.02,0.02))+
   scale_color_manual(values=my_colours) +
   scale_fill_manual(values=my_colours)
   Model4b_plot

  #ggsave(file="Model4b_plot.png", plot=Model4b_plot, width=4, height=4)
  #ggsave(file="Model4b_plot.pdf", plot=Model4b_plot, width=4, height=4)

Model 5b: Time-interval 4 - ENSO.98_14 with lme4

  # 1. Build lme4 model:
   
    model5b <- lmer(
      Coral_cover ~ Thermal_regime * Year + 
        (1|Region), data=ENSO.83_14, REML = FALSE)
  
    summary(model5b)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: Coral_cover ~ Thermal_regime * Year + (1 | Region)
##    Data: ENSO.83_14
## 
##      AIC      BIC   logLik deviance df.resid 
##   1443.8   1475.1   -711.9   1423.8      159 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.52432 -0.48425 -0.05873  0.47939  3.06957 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Region   (Intercept)  94.08    9.699  
##  Residual             234.31   15.307  
## Number of obs: 169, groups:  Region, 13
## 
## Fixed effects:
##                                           Estimate Std. Error         df
## (Intercept)                             -1737.0757   506.0616   162.8735
## Thermal_regimeTropical upwelling         1854.0315   708.3227   167.0809
## Thermal_regimeEquatorial upwelling        291.8845  1024.3367   165.7040
## Thermal_regimeSeasonal                   2738.6990   918.4800   167.5534
## Year                                        0.8799     0.2530   162.6103
## Thermal_regimeTropical upwelling:Year      -0.9207     0.3539   166.9356
## Thermal_regimeEquatorial upwelling:Year    -0.1475     0.5117   165.5744
## Thermal_regimeSeasonal:Year                -1.3650     0.4588   167.3530
##                                         t value Pr(>|t|)    
## (Intercept)                              -3.433 0.000758 ***
## Thermal_regimeTropical upwelling          2.617 0.009671 ** 
## Thermal_regimeEquatorial upwelling        0.285 0.776038    
## Thermal_regimeSeasonal                    2.982 0.003293 ** 
## Year                                      3.478 0.000648 ***
## Thermal_regimeTropical upwelling:Year    -2.601 0.010116 *  
## Thermal_regimeEquatorial upwelling:Year  -0.288 0.773525    
## Thermal_regimeSeasonal:Year              -2.975 0.003365 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Thr_Tu Thr_Eu Thrm_S Year   T_Tu:Y T_Eu:Y
## Thrml_rgmTu -0.714                                          
## Thrml_rgmEu -0.494  0.353                                   
## Thrml_rgmSs -0.537  0.383  0.265                            
## Year        -1.000  0.714  0.494  0.536                     
## Thrml_rTu:Y  0.715 -1.000 -0.353 -0.383 -0.715              
## Thrml_rEu:Y  0.494 -0.353 -1.000 -0.265 -0.494  0.353       
## Thrml_rgS:Y  0.537 -0.383 -0.265 -1.000 -0.536  0.383  0.265
    anova(model5b) # Thermal regime and its interaciton with year are significant 
## Type III Analysis of Variance Table with Satterthwaite's method
##                      Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## Thermal_regime      2922.20  974.07     3 167.16  4.1572 0.007168 **
## Year                 580.02  580.02     1 166.90  2.4755 0.117530   
## Thermal_regime:Year 2898.06  966.02     3 167.02  4.1228 0.007496 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ranova(model5b) # ramdom effects are significant
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## Coral_cover ~ Thermal_regime + Year + (1 | Region) + Thermal_regime:Year
##              npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>         10 -711.89 1443.8                         
## (1 | Region)    9 -718.96 1455.9 14.146  1  0.0001692 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    drop1(model5b)
## Single term deletions using Satterthwaite's method:
## 
## Model:
## Coral_cover ~ Thermal_regime * Year + (1 | Region)
##                     Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## Thermal_regime:Year 2898.1  966.02     3 167.02  4.1228 0.007496 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  # 2. Predict values:
    pred5 <- predict(model5b,re.form = NA)

  #3. Bootstrap CI:
    boot5 <- bootMer(model5b, predict, nsim = 1000, re.form = NULL) # Include ramndom effects
    std.err <- apply(boot5$t, 2, sd)
    CI.lo_5 <- pred5 - std.err*1.96
    CI.hi_5 <- pred5 + std.err*1.96

  #Plot
   Model5b_plot<- ggplot(ENSO.83_14, aes(
     x = Year, y = Coral_cover, colour = Thermal_regime)) +
    geom_line(aes(y = pred5),size=2) +
    geom_point(aes(fill=factor(Thermal_regime)),
             shape = 21, colour = "black", size = 2, 
             stroke = 0.3, alpha=0.5) +
    geom_ribbon(aes(ymin = CI.lo_5, ymax = CI.hi_5),size=2, 
                alpha = 0.03, linetype = 0) +
    scale_y_continuous("Live coral cover (%)",
                       limits = c(-20, 90), 
                       breaks = seq(0, 90, by=20), expand = c(0,0))+
    scale_x_continuous("", limits = c(1983, 2014.3),
                     breaks = seq(1983, 2014, by=2), expand = c(0.02,0.02))+
   scale_color_manual(values=my_colours) +
   scale_fill_manual(values=my_colours)
   Model5b_plot

  #ggsave(file="Model4b_plot.png", plot=Model4b_plot, width=4, height=4)
  #ggsave(file="Model4b_plot.pdf", plot=Model4b_plot, width=4, height=4)